Supervised Learning
Unsupervised Learning
Secara umum, unsupervised learning yang akan kita pelajari terbagi menjadi dua:
Tujuan dimensionality reduction adalah untuk mereduksi banyaknya variabel (dimensi/fitur) pada data dengan tetap mempertahankan informasi sebanyak mungkin. Dimensionality reduction dapat mengatasi masalah high-dimensional data. Kesulitan yang dihadapi pada high-dimensional data:
Istilah:
💡 Mari kita ulas sedikit terkait materi Practical Statistics yaitu variance:
a <- c(1,2,3,4,5,6,7,2,9,10)
b <- c(4,5,5,6,6,4,6,5,4,4)
stripchart(a, method = "jitter")stripchart(b, method = "jitter")var(a)#> [1] 9.433333
var(b)#> [1] 0.7666667
Tapi, hati-hati nilai variance bergantung pada range dari variabel.
x <- a*100
x#> [1] 100 200 300 400 500 600 700 200 900 1000
var(x)#> [1] 94333.33
Notes:
Informasi adalah nilai variance yang sudah discaling
Pada data gambar, setiap kotak pixel akan menjadi 1 kolom. Foto berukuran 40x40 pixel memiliki 1600 kolom (dimensi). Sekarang mari renungkan, berapa spesifikasi kamera handphone anda? Berapa besar dimensi data yang dihasilkan kamera Anda?
Image compression adalah salah satu penerapan dimensionality reduction. Data gambar direduksi dimensinya namun tetap menghasilkan gambar yang serupa (informasi inti tidak hilang), sehingga gambar lebih mudah diproses. Salah satu algoritma dimensionality reduction adalah Principal Component Analysis (PCA).
Untuk lebih detail terkait image compression, silahkan kunjungi artikel dari salah satu tim kami: https://medium.com/analytics-vidhya/image-compression-using-k-means-clustering-and-principal-component-analysis-in-python-893fd37e8e15
Ide dasar dari PCA adalah untuk membuat sumbu (axis) baru yang dapat menangkap informasi sebesar mungkin. Sumbu baru ini adalah yang dinamakan sebagai Principal Component (PC). Untuk melakukan dimensionality reduction, kita akan memilih beberapa PC untuk dapat merangkum informasi yang dibutuhkan.
Figure 1A (kiri):
Figure 1B (kanan):
Hanya dengan menggunakan 1 PC (yaitu PC1), kita dapat mereduksi dimensi data sebesar 50% (dari 2 dimensi menjadi 1 dimensi), namun tetap mempertahankan informasi sebesar 90% informasi data asli.
💡 Catatan:
Di antara kedua data berikut, mana yang lebih cocok untuk dilakukan PCA:
Tujuan PCA: untuk reduksi dimensi, kira-kira data yang seperti apa yang tepat digunakan untuk reduksi dimensi menggunakan PCA?
Jawaban: Logistic Machinery
💡 Kunjungi link berikut untuk ilustrasi secara visual untuk data 2D dan 3D: https://setosa.io/ev/principal-component-analysis/
⚠️ Miskonsepsi garis PCA dengan garis Linear Regression
Berikut adalah animasi bagaimana garis PCA secara intuitif dibentuk:
Source: Medium
Notes:
PCA bertujuan untuk dimensionality reduction dengan cara fitur extraction (mengekstrak fitur yang ada)
x1,x2,x3 —> PC1, PC2, PC3
x1, x2, x3,x4, … , x10 —-> PC1, PC2, … , PC10
yang mana PC1 (variance paling banyak) 70% PC2 (variance menengah jumlahnya) 20% PC3 (variance paling kecil) 10% PC4 (variance lebih kecil dr PC3) 0%
PC1 = 0.4x1 + 0.3x2 + 0.2 x3 + 0.1 x4 PC2 = …x1 + ..x2 + …
Misal kita punya 4 variabel (X1, X2, X3, X4), kemudian kita lakukan PCA sehingga kita memiliki 4 variabel baru (PC1, PC2, PC3, PC4). Ketika kita ingin mempertahankan minimal 80% informasi, berapa PC/dimensi yang kita kurangi?
PC1 (variance paling banyak) 70% PC2 (variance menengah jumlahnya) 20% PC3 (variance paling kecil) 10% PC4 (variance lebih kecil dr PC3) 0%
Jadi PCA akan mengubah data asli menjadi data proyeksinya
Kita sebagai data scientist diminta untuk melakukan dimensionality reduction terhadap data unit bangunan yang terjual di New York Property Market selama periode 12 bulan. Sumber data diperoleh dari New York City Department of Finance.
library(dplyr)
property <- read.csv("data_input/nyc.csv")
glimpse(property)#> Rows: 84,548
#> Columns: 22
#> $ X <int> 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1~
#> $ BOROUGH <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ~
#> $ NEIGHBORHOOD <chr> "ALPHABET CITY", "ALPHABET CITY", "ALPH~
#> $ BUILDING.CLASS.CATEGORY <chr> "07 RENTALS - WALKUP APARTMENTS ~
#> $ TAX.CLASS.AT.PRESENT <chr> "2A", "2", "2", "2B", "2A", "2", "2B", ~
#> $ BLOCK <int> 392, 399, 399, 402, 404, 405, 406, 407,~
#> $ LOT <int> 6, 26, 39, 21, 55, 16, 32, 18, 34, 153,~
#> $ EASE.MENT <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,~
#> $ BUILDING.CLASS.AT.PRESENT <chr> "C2", "C7", "C7", "C4", "C2", "C4", "C4~
#> $ ADDRESS <chr> "153 AVENUE B", "234 EAST 4TH STREET"~
#> $ APARTMENT.NUMBER <chr> " ", " ", " ", " ", " ", " ", " ", " ",~
#> $ ZIP.CODE <int> 10009, 10009, 10009, 10009, 10009, 1000~
#> $ RESIDENTIAL.UNITS <int> 5, 28, 16, 10, 6, 20, 8, 44, 15, 24, 30~
#> $ COMMERCIAL.UNITS <int> 0, 3, 1, 0, 0, 0, 0, 2, 0, 0, 4, 0, 0, ~
#> $ TOTAL.UNITS <int> 5, 31, 17, 10, 6, 20, 8, 46, 15, 24, 34~
#> $ LAND.SQUARE.FEET <chr> "1633", "4616", "2212", "2272", "2369",~
#> $ GROSS.SQUARE.FEET <chr> "6440", "18690", "7803", "6794", "4615"~
#> $ YEAR.BUILT <int> 1900, 1900, 1900, 1913, 1900, 1900, 192~
#> $ TAX.CLASS.AT.TIME.OF.SALE <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, ~
#> $ BUILDING.CLASS.AT.TIME.OF.SALE <chr> "C2", "C7", "C7", "C4", "C2", "C4", "C4~
#> $ SALE.PRICE <chr> "6625000", " - ", " - ", "3936272", "~
#> $ SALE.DATE <chr> "2017-07-19 00:00:00", "2016-12-14 00:0~
Deskripsi data:
X: ID barisBOROUGH: Kode digit untuk wilayah dimana properti berada; Manhattan (1), Bronx (2), Brooklyn (3), Queens (4), dan Staten Island (5)NEIGHBORHOOD: Nama lingkunganBUILDING.CLASS.CATEGORY: Kategori kelas propertiTAX.CLASS.AT.PRESENT, TAX.CLASS.AT.TIME.OF.SALE: Jenis pajak bangunanBLOCK, LOT: Kombinasi borough, blok, dan lot membentuk kode unik untuk properti di New York CityEASE.MENT: Ragam hak yang dimiliki bangunan (contoh: hak jalan, dsb.)BUILDING.CLASS.AT.PRESENT, BUILDING.CLASS.AT.TIME.OF.SALE: Kategori kepemilikan bangunan (A: rumah satu keluarga; O: gedung kantor; R: kondominium (kepemilikan bersama); dsb.)ADDRESS: Alamat propertiAPARTMENT.NUMBER: Nomor apartemenZIP.CODE: Kode pos propertiRESIDENTIAL.UNITS: Jumlah unit hunian di properti yang terdaftarCOMMERCIAL.UNITS: Jumlah unit komersial di properti yang terdaftarTOTAL.UNITS: Jumlah unit hunian dan komersialLAND.SQUARE.FEET: Luas tanah properti (square feet)GROSS.SQUARE.FEET: Total area semua lantai bangunan yang diukur dari permukaan luar dinding luar bangunan, termasuk area tanah dan ruang di dalam setiap bangunan atau struktur pada propertiYEAR.BUILT: Tahun properti dibangunSALE.PRICE: Harga properti terjual; bila 0 maka properti merupakan warisanSALE.DATE: Tanggal properti dijualBuang kolom yang tidak diperlukan:
X (hanya ID baris)EASE.MENT (semua NA)Memperbaiki tipe data yang belum tepat:
Perhatikan yang integer:
BOROUGH, BLOCK, LOT: int -> factorZIP.CODE: int -> factorRESIDENTIAL.UNITS, COMMERCIAL.UNITS, TOTAL.UNITS: int -> int (sudah benar)YEAR.BUILT: int -> int (tapi bisa faktor, tergantung analisa yang digunakan)TAX.CLASS.AT.TIME.OF.SALE: int -> factorPerhatikan yang character:
NEIGHBORHOOD: chr -> fctBUILDING.CLASS.CATEGORY, BUILDING.CLASS.AT.PRESENT, BUILDING.CLASS.AT.TIME.OF.SALE: chr -> fctTAX.CLASS.AT.PRESENT: chr -> fctADDRESS: chr -> chr (sudah benar)APARTMENT.NUMBER: chr -> chr (sudah benar)LAND.SQUARE.FEET, GROSS.SQUARE.FEET: chr -> numSALE.PRICE: chr -> numSALE.DATE: chr -> date/ posixct (date and time)library(lubridate)
property_clean <- property %>%
select(-c(X, EASE.MENT)) %>%
mutate_at(.vars = c("BOROUGH", "BLOCK", "LOT", "ZIP.CODE", "TAX.CLASS.AT.TIME.OF.SALE", "NEIGHBORHOOD", "BUILDING.CLASS.CATEGORY", "BUILDING.CLASS.AT.PRESENT", 'BUILDING.CLASS.AT.TIME.OF.SALE', "TAX.CLASS.AT.PRESENT"), as.factor) %>%
mutate_at(.vars = c("LAND.SQUARE.FEET", "GROSS.SQUARE.FEET", "SALE.PRICE"),
as.numeric) %>%
mutate(SALE.DATE = ymd_hms(SALE.DATE))
head(property_clean)Cek missing values:
property_clean %>%
is.na() %>%
colSums()#> BOROUGH NEIGHBORHOOD
#> 0 0
#> BUILDING.CLASS.CATEGORY TAX.CLASS.AT.PRESENT
#> 0 0
#> BLOCK LOT
#> 0 0
#> BUILDING.CLASS.AT.PRESENT ADDRESS
#> 0 0
#> APARTMENT.NUMBER ZIP.CODE
#> 0 0
#> RESIDENTIAL.UNITS COMMERCIAL.UNITS
#> 0 0
#> TOTAL.UNITS LAND.SQUARE.FEET
#> 0 26252
#> GROSS.SQUARE.FEET YEAR.BUILT
#> 27612 0
#> TAX.CLASS.AT.TIME.OF.SALE BUILDING.CLASS.AT.TIME.OF.SALE
#> 0 0
#> SALE.PRICE SALE.DATE
#> 14561 0
Kolom dengan missing values: GROSS.SQUARE.FEET, LAND.SQUARE.FEET , SALE.PRICE
Beberapa cara untuk handle missing values:
Catatan: Fokus pada course ini bukan di handle NA, namun untuk memahami PCA dan visualisasinya, sehingga mari kita gunakan cara yang cukup simpel, yaitu cara ketiga. Dengan justifikasi bahwa kita tetap membutuhkan informasi dari kolom tersebut dalam analisis kita.
# drop baris yang mengandung NA
property_clean <- property_clean %>%
filter(complete.cases(.)) # complete cases utk filter data-data yang barisnya tidak ada NA
anyNA(property_clean)#> [1] FALSE
dim(property_clean)#> [1] 48244 20
head(property_clean)glimpse(property_clean)#> Rows: 48,244
#> Columns: 20
#> $ BOROUGH <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ~
#> $ NEIGHBORHOOD <fct> ALPHABET CITY, ALPHABET CITY, ALPHABET ~
#> $ BUILDING.CLASS.CATEGORY <fct> 07 RENTALS - WALKUP APARTMENTS ~
#> $ TAX.CLASS.AT.PRESENT <fct> 2A, 2B, 2A, 2B, 2, 2B, 2, 2A, 2A, 2A, 4~
#> $ BLOCK <fct> 392, 402, 404, 406, 387, 400, 376, 391,~
#> $ LOT <fct> 6, 21, 55, 32, 153, 21, 14, 19, 4, 5, 3~
#> $ BUILDING.CLASS.AT.PRESENT <fct> C2, C4, C2, C4, D9, D1, C6, S3, S4, S5,~
#> $ ADDRESS <chr> "153 AVENUE B", "154 EAST 7TH STREET", ~
#> $ APARTMENT.NUMBER <chr> " ", " ", " ", " ", " ", " ", " ", " ",~
#> $ ZIP.CODE <fct> 10009, 10009, 10009, 10009, 10009, 1000~
#> $ RESIDENTIAL.UNITS <int> 5, 10, 6, 8, 24, 10, 24, 3, 4, 5, 0, 1,~
#> $ COMMERCIAL.UNITS <int> 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, ~
#> $ TOTAL.UNITS <int> 5, 10, 6, 8, 24, 10, 24, 4, 5, 6, 1, 1,~
#> $ LAND.SQUARE.FEET <dbl> 1633, 2272, 2369, 1750, 4489, 3717, 413~
#> $ GROSS.SQUARE.FEET <dbl> 6440, 6794, 4615, 4226, 18523, 12350, 1~
#> $ YEAR.BUILT <int> 1900, 1913, 1900, 1920, 1920, 2009, 192~
#> $ TAX.CLASS.AT.TIME.OF.SALE <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 1, 1, ~
#> $ BUILDING.CLASS.AT.TIME.OF.SALE <fct> C2, C4, C2, C4, D9, D1, C6, S3, S4, S5,~
#> $ SALE.PRICE <dbl> 6625000, 3936272, 8000000, 3192840, 162~
#> $ SALE.DATE <dttm> 2017-07-19, 2016-09-23, 2016-11-17, 20~
Dikarenakan analisis PCA menggunakan nilai variance, kita hanya dapat mereduksi kolom bertipe data numerik.
# data untuk PCA
ppt <- property_clean %>%
select_if(is.numeric)
head(ppt)Apakah antar variable memiliki skala yang sama? Mari kita cek range untuk tiap kolom:
summary(ppt)#> RESIDENTIAL.UNITS COMMERCIAL.UNITS TOTAL.UNITS LAND.SQUARE.FEET
#> Min. : 0.000 Min. : 0.0000 Min. : 0.000 Min. : 0
#> 1st Qu.: 1.000 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.: 1413
#> Median : 1.000 Median : 0.0000 Median : 1.000 Median : 2140
#> Mean : 2.566 Mean : 0.2492 Mean : 2.835 Mean : 3358
#> 3rd Qu.: 2.000 3rd Qu.: 0.0000 3rd Qu.: 2.000 3rd Qu.: 3071
#> Max. :1844.000 Max. :2261.0000 Max. :2261.000 Max. :4228300
#> GROSS.SQUARE.FEET YEAR.BUILT SALE.PRICE
#> Min. : 0 Min. : 0 Min. : 0
#> 1st Qu.: 828 1st Qu.:1920 1st Qu.: 80420
#> Median : 1620 Median :1931 Median : 480000
#> Mean : 3670 Mean :1828 Mean : 1153281
#> 3rd Qu.: 2520 3rd Qu.:1961 3rd Qu.: 830000
#> Max. :3750565 Max. :2017 Max. :2210000000
Kesimpulan: beda
Mari kita cek juga matriks covariance dari ppt:
cov(ppt)#> RESIDENTIAL.UNITS COMMERCIAL.UNITS TOTAL.UNITS
#> RESIDENTIAL.UNITS 305.043041 2.451582 307.4394
#> COMMERCIAL.UNITS 2.451582 120.737718 123.1811
#> TOTAL.UNITS 307.439418 123.181081 430.5798
#> LAND.SQUARE.FEET 220192.076495 18262.941834 238387.9278
#> GROSS.SQUARE.FEET 318357.747484 22509.684842 340793.2306
#> YEAR.BUILT 220.658236 27.776918 242.6401
#> SALE.PRICE 28718401.580456 6953367.278400 35653380.3633
#> LAND.SQUARE.FEET GROSS.SQUARE.FEET YEAR.BUILT
#> RESIDENTIAL.UNITS 220192.08 318357.75 220.65824
#> COMMERCIAL.UNITS 18262.94 22509.68 27.77692
#> TOTAL.UNITS 238387.93 340793.23 242.64015
#> LAND.SQUARE.FEET 988215601.32 547933690.70 138854.30577
#> GROSS.SQUARE.FEET 547933690.70 868770428.70 429818.63540
#> YEAR.BUILT 138854.31 429818.64 215631.28007
#> SALE.PRICE 17537208440.02 179859907831.20 47709930.20464
#> SALE.PRICE
#> RESIDENTIAL.UNITS 28718402
#> COMMERCIAL.UNITS 6953367
#> TOTAL.UNITS 35653380
#> LAND.SQUARE.FEET 17537208440
#> GROSS.SQUARE.FEET 179859907831
#> YEAR.BUILT 47709930
#> SALE.PRICE 179595088101951
Variance dari masing-masing variabel berbeda jauh karena range/skala dari tiap variabel berbeda, begitupun covariance. Nilai variance dan covariance dipengaruhi oleh skala dari data. Semakin tinggi skala, nilai variance atau covariance akan semakin tinggi.
⚠️ Data dengan perbedaan skala antar variabel yang tinggi tidak baik untuk langsung dianalisis PCA karena dapat menimbulkan bias. PC1 akan dianggap menangkap variansi tertinggi dan PC selanjutnya dianggap tidak memberikan informasi.
Perhatikan plot di bawah ini bahwa semua informasi hanya dirangkum oleh PC1. Hal ini terjadi karena SALE.PRICE memiliki skala hingga milyaran (variansinya tinggi) dibandingkan variabel lainnya.
💡 Oleh sebab itu, harus kita scaling terlebih dahulu sebelum melakukan PCA.
Scaling dilakukan agar antar variabel memiliki skala yang tidak jauh berbeda. Nilai scaling yang digunakan adalah Z-score (mean = 0, standar deviasi = 1)
\[Z = \frac{x-mean}{sd}\]
# scaling
ppt_z <- scale(ppt)
head(ppt)head(ppt_z)#> RESIDENTIAL.UNITS COMMERCIAL.UNITS TOTAL.UNITS LAND.SQUARE.FEET
#> [1,] 0.1393299 -0.02267649 0.1043530 -0.05487728
#> [2,] 0.4256088 -0.02267649 0.3453120 -0.03455020
#> [3,] 0.1965857 -0.02267649 0.1525448 -0.03146456
#> [4,] 0.3110972 -0.02267649 0.2489284 -0.05115542
#> [5,] 1.2271899 -0.02267649 1.0199972 0.03597427
#> [6,] 0.4256088 -0.02267649 0.3453120 0.01141636
#> GROSS.SQUARE.FEET YEAR.BUILT SALE.PRICE
#> [1,] 0.09398662 0.1555574 0.4082973
#> [2,] 0.10599684 0.1835529 0.2076656
#> [3,] 0.03206956 0.1555574 0.5108992
#> [4,] 0.01887190 0.1986274 0.1521910
#> [5,] 0.50392848 0.1986274 1.1251676
#> [6,] 0.29449613 0.3902885 0.6862552
# melihat variance yang dirangkum tiap PC
plot(prcomp(ppt_z))Sekarang, PC1 sudah tidak mendominasi PC lainnya dalam besarnya variance yang ditangkap.
prcomp()PCA dapat dilakukan dengan menggunakan fungsi prcomp().
Cara 1: menggunakan data yang di-scale secara terpisah, yaitu ppt_z
# lakukan scaling (tadi sudah dilakukan, disimpan dalam object ppt_z)
# buat PCA dgn prcomp()
prcomp(ppt_z)#> Standard deviations (1, .., p=7):
#> [1] 1.691768718 1.094108600 1.007258872 0.999513415 0.842647844 0.466010414
#> [7] 0.005136076
#>
#> Rotation (n x k) = (7 x 7):
#> PC1 PC2 PC3 PC4 PC5
#> RESIDENTIAL.UNITS -0.50181439 0.03289418 0.265231639 -0.001822824 0.52716521
#> COMMERCIAL.UNITS -0.18992921 -0.75302340 -0.358571533 -0.003172376 -0.34847875
#> TOTAL.UNITS -0.52285631 -0.37109583 0.033364411 -0.003808758 0.25919705
#> LAND.SQUARE.FEET -0.38007544 0.21543634 0.357103052 -0.038749370 -0.71562101
#> GROSS.SQUARE.FEET -0.50150909 0.31906907 -0.102415876 -0.011363727 -0.13797633
#> YEAR.BUILT -0.02620196 0.01241830 0.003251844 0.999110532 -0.02805113
#> SALE.PRICE -0.20512249 0.38182393 -0.813598964 -0.010936303 0.03922413
#> PC6 PC7
#> RESIDENTIAL.UNITS 0.20660513 -0.59679168068
#> COMMERCIAL.UNITS -0.07659666 -0.37549241058
#> TOTAL.UNITS 0.13324195 0.70911561460
#> LAND.SQUARE.FEET 0.40986865 0.00005036682
#> GROSS.SQUARE.FEET -0.78551126 -0.00001709229
#> YEAR.BUILT 0.01182229 0.00042438145
#> SALE.PRICE 0.38540088 0.00004079992
Cara 2: menggunakan data yang belum discale, lalu tambahkan parameter scale = T
# memasukkan data yang belum discaling ke dalam fungsi prcomp, dan tambahkan scale =T
pca <- prcomp(ppt, scale = T)Terdapat tiga komponen utama dalam objek pca:
pca$sdev: standar deviasi (akar variance) yang ditangkap oleh masing-masing PC, disebut juga sebagai eigen value. Digunakan untuk mengetahui seberapa besaran informasi yang ditangkap masing-masing PC.pca$sdev # akar dr variance#> [1] 1.691768718 1.094108600 1.007258872 0.999513415 0.842647844 0.466010414
#> [7] 0.005136076
# variance yang ditangkap, disebut juga eigen value
pca$sdev^2 #> [1] 2.86208139528 1.19707362796 1.01457043534 0.99902706733 0.71005538886
#> [6] 0.21716570595 0.00002637928
PC1 akan secara otomatis merangkum paling banyak informasi (nilai variancenya paling besar), kemudian diikuti oleh PC2, PC3, …, PC7
pca$rotation: matrix rotasi yang berfungsi untuk memproyeksikan titik ke masing-masing PC, matrix ini terdiri dari eigen vector. Digunakan untuk mengetahui kontribusi masing-masing variabel ke PC.pca$rotation#> PC1 PC2 PC3 PC4 PC5
#> RESIDENTIAL.UNITS -0.50181439 0.03289418 0.265231639 -0.001822824 0.52716521
#> COMMERCIAL.UNITS -0.18992921 -0.75302340 -0.358571533 -0.003172376 -0.34847875
#> TOTAL.UNITS -0.52285631 -0.37109583 0.033364411 -0.003808758 0.25919705
#> LAND.SQUARE.FEET -0.38007544 0.21543634 0.357103052 -0.038749370 -0.71562101
#> GROSS.SQUARE.FEET -0.50150909 0.31906907 -0.102415876 -0.011363727 -0.13797633
#> YEAR.BUILT -0.02620196 0.01241830 0.003251844 0.999110532 -0.02805113
#> SALE.PRICE -0.20512249 0.38182393 -0.813598964 -0.010936303 0.03922413
#> PC6 PC7
#> RESIDENTIAL.UNITS 0.20660513 -0.59679168068
#> COMMERCIAL.UNITS -0.07659666 -0.37549241058
#> TOTAL.UNITS 0.13324195 0.70911561460
#> LAND.SQUARE.FEET 0.40986865 0.00005036682
#> GROSS.SQUARE.FEET -0.78551126 -0.00001709229
#> YEAR.BUILT 0.01182229 0.00042438145
#> SALE.PRICE 0.38540088 0.00004079992
PC1 = -0.50181439 x residential.unit - 0.18992921 x commercial.unit +…
baris 1 di data awal = 2 residential unit, 3 commercial unit + … baris 1 untuk PC1 = -0.50181439 x 2 - 0.18992921 x 3 +… = xxx
pca$x: nilai hasil proyeksi titik ke PC untuk tiap baris. Digunakan untuk mendapatkan nilai data yang baru.as.data.frame(ppt_z)as.data.frame(pca$x)Review Day 1:
Dimensionality reduction: - Feature selection = menseleksi fitur/variable dengan cara membuang variabel2 yang sekiranya tidak berpengaruh (misal: awal 10 variabel –> kita hapus 2 var –> sisa 8 var). Contoh metode: menghapus var yg berkorelasi sgt tinggi, nearzerovar (membuang var yg variansinya ~0) - Feature extraction = mengekstrak fitur/ variabel. Caranya adalah dengan membobotkan semua variabel ke dalam variabel baru (misal: awal 10 var –> 10 var baru dengan value yg berbeda –> dibuang beberapa var yang variansinya rendah/ memilih PC). metode: PCA
Fungsi untuk membuat PCA (tanpa library):
pca <- prcomp(data_numerik, scale = T)
pca$sdev^2 # variance adl kuadrat sdev (variance = eigen value)
pca$rotation
pca$x
Workflow:
head(property_clean)Pilih var yang numeric saja
ppt <- property_clean %>%
select_if(is.numeric)Cek skala data, apakah sudah sama?
summary(ppt)#> RESIDENTIAL.UNITS COMMERCIAL.UNITS TOTAL.UNITS LAND.SQUARE.FEET
#> Min. : 0.000 Min. : 0.0000 Min. : 0.000 Min. : 0
#> 1st Qu.: 1.000 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.: 1413
#> Median : 1.000 Median : 0.0000 Median : 1.000 Median : 2140
#> Mean : 2.566 Mean : 0.2492 Mean : 2.835 Mean : 3358
#> 3rd Qu.: 2.000 3rd Qu.: 0.0000 3rd Qu.: 2.000 3rd Qu.: 3071
#> Max. :1844.000 Max. :2261.0000 Max. :2261.000 Max. :4228300
#> GROSS.SQUARE.FEET YEAR.BUILT SALE.PRICE
#> Min. : 0 Min. : 0 Min. : 0
#> 1st Qu.: 828 1st Qu.:1920 1st Qu.: 80420
#> Median : 1620 Median :1931 Median : 480000
#> Mean : 3670 Mean :1828 Mean : 1153281
#> 3rd Qu.: 2520 3rd Qu.:1961 3rd Qu.: 830000
#> Max. :3750565 Max. :2017 Max. :2210000000
Kalau belum sama, lakukan scaling
pca <- prcomp(ppt, scale = T)Cek variansinya/ eigen value
pca$sdev^2#> [1] 2.86208139528 1.19707362796 1.01457043534 0.99902706733 0.71005538886
#> [6] 0.21716570595 0.00002637928
cek bobot/ kontribusi tiap var thd PC
pca$rotation#> PC1 PC2 PC3 PC4 PC5
#> RESIDENTIAL.UNITS -0.50181439 0.03289418 0.265231639 -0.001822824 0.52716521
#> COMMERCIAL.UNITS -0.18992921 -0.75302340 -0.358571533 -0.003172376 -0.34847875
#> TOTAL.UNITS -0.52285631 -0.37109583 0.033364411 -0.003808758 0.25919705
#> LAND.SQUARE.FEET -0.38007544 0.21543634 0.357103052 -0.038749370 -0.71562101
#> GROSS.SQUARE.FEET -0.50150909 0.31906907 -0.102415876 -0.011363727 -0.13797633
#> YEAR.BUILT -0.02620196 0.01241830 0.003251844 0.999110532 -0.02805113
#> SALE.PRICE -0.20512249 0.38182393 -0.813598964 -0.010936303 0.03922413
#> PC6 PC7
#> RESIDENTIAL.UNITS 0.20660513 -0.59679168068
#> COMMERCIAL.UNITS -0.07659666 -0.37549241058
#> TOTAL.UNITS 0.13324195 0.70911561460
#> LAND.SQUARE.FEET 0.40986865 0.00005036682
#> GROSS.SQUARE.FEET -0.78551126 -0.00001709229
#> YEAR.BUILT 0.01182229 0.00042438145
#> SALE.PRICE 0.38540088 0.00004079992
baris 1 PC1 = -0.50181439 x baris 1 residential unit -0.18992921 x baris 1 commercial unit + … baris 2 PC1 =
cek data baru hasil PCA
pca$x %>%
as.data.frame()Misal kita hanya ingin mereduksi 7 dimensi menjadi 5 dimensi atau saya ingin mempertahankan minimal 90% informasi, maka kita gunakan kolom PC1 sampai PC5 pada pca$x. Namun bagaimana kita memilih berapa dimensi yang kita pertahankan? Kita melihat cumulative variance untuk melakukan dimensionality reduction dengan fungsi summary().
summary(pca)#> Importance of components:
#> PC1 PC2 PC3 PC4 PC5 PC6 PC7
#> Standard deviation 1.6918 1.0941 1.0073 0.9995 0.8426 0.46601 0.005136
#> Proportion of Variance 0.4089 0.1710 0.1449 0.1427 0.1014 0.03102 0.000000
#> Cumulative Proportion 0.4089 0.5799 0.7248 0.8675 0.9690 1.00000 1.000000
pca$sdevPemilihan banyaknya PC disesuaikan dengan kebutuhan informasi. Misal, bila kita ingin merangkum 75% informasi, maka jumlah PC yang kita gunakan adalah 4 (PC1 sampai PC4)
pc_keep <- pca$x[,1:4] %>%
as.data.frame()
pc_keepCek perbedaan sebelum dilakukan PCA dan setelah
ppt_z %>%
as.data.frame()pca$x %>%
as.data.frame()Setelah dipilih PC yang merangkum informasi yang dibutuhkan, PC dapat digabung dengan data awal dan digunakan untuk analisis lebih lanjut (misal: supervised learning).
property_clean %>%
select_if(purrr::negate(is.numeric)) %>% # ambil kolom selain numeric
cbind(pc_keep) # gabungkan dengan PC1 - PC4Kelebihan melakukan PCA sebelum modeling (regresi/klasifikasi/clustering):
library(GGally)
# korelasi sebelum PCA
ggcorr(ppt, label = T, hjust = 1)# korelasi setelah PCA
ggcorr(pc_keep, label = T)Additional: kalau var sale price ingin digunakan sebagai target variabel agar gak masuk ke PCA
ppt <- property_clean %>%
select(-SALE.PRICE) %>%
select_if(is.numeric)Kekurangan melakukan PCA sebelum modeling:
💡 Video: 3Blue1Brown: Eigenvectors and eigenvalues
Untuk membentuk PC dan nilai pada PC dibutuhkan eigen values & eigen vector. Secara manual, eigen values dan eigen vector didapatkan dari operasi matrix.
Teori matrix:
Perkalian skalar-vektor: mengubah besaran vektor (tidak merubah arah kecuali berbalik arah)
\[2\left(\begin{array}{cc}2\\3\end{array}\right) = \left(\begin{array}{cc}4\\6\end{array}\right)\] \[-2\left(\begin{array}{cc}2\\3\end{array}\right) = \left(\begin{array}{cc}-4\\-6\end{array}\right)\]
Perkalian matrix-vektor: mengubah besaran dan arah
\[\left(\begin{array}{cc}1 & 2\\-1 & 2\end{array}\right)\left(\begin{array}{cc}2\\3 \end{array}\right) = \left(\begin{array}{cc}8\\4\end{array}\right)\]
Perkalian Matrix-Vektor Spesial
\[ \left(\begin{array}{cc} 0 & -1\\ 1 & 0 \end{array}\right) \left(\begin{array}{cc} 2\\ 3 \end{array}\right)= \left(\begin{array}{cc} -3\\ 2 \end{array}\right) \]
# perkalian matrix di R
matrix(c(0,1,-1,0), nrow=2) %*% as.vector(c(2,3))#> [,1]
#> [1,] -3
#> [2,] 2
\[ \left(\begin{array}{cc} 1 & 0\\ 0 & 1 \end{array}\right) \left(\begin{array}{cc} 2\\ 3 \end{array}\right)= \left(\begin{array}{cc} 2\\ 3 \end{array}\right) \]
# perkalian matrix di R
matrix(c(1,0,0,1), nrow=2) %*% as.vector(c(2,3))#> [,1]
#> [1,] 2
#> [2,] 3
Eigen dari suatu Matrix
Untuk setiap matrix, terdapat vektor spesial (eigen vector) yang jika dikalikan dengan matrixnya, hasilnya akan sama dengan vektor tersebut dikalikan suatu skalar (eigen value). Sehingga didapatkan rumus:
\[Ax = \lambda x\] Dimana \(x\) adalah eigen vector dan \(\lambda\) adalah eigen value dari matrix \(A\).
\[ \left(\begin{array}{cc} 2 & 3\\ 2 & 1 \end{array}\right) \left(\begin{array}{cc} 3\\ 2 \end{array}\right) = \left(\begin{array}{cc} 12\\ 8 \end{array}\right) =4 \left(\begin{array}{cc} 3\\ 2 \end{array}\right) \]
Teori eigen dipakai untuk menentukan PC dan nilai-nilai pada PC.
Penerapan Eigen dalam PCA:
Matrix covariance adalah matrix yang dapat merangkum informasi (variance) dari data. Kita menggunakan matrix covariance untuk mendapatkan mendapatkan eigen vector dan eigen value dari matrix tersebut, dimana:
Eigen vector akan menjadi formula untuk kalkulasi nilai di setiap PC. Contohnya, untuk data yang terdiri dari 2 variabel, bila diketahui eigen vector dari PC1 adalah:
\[x_{PC1}= \left(\begin{array}{cc}b_1\\b_2\end{array}\right)\]
Maka formula untuk menghitung nilai pada PC1 (untuk tiap barisnya) adalah:
\[PC1= b_1X1 + b_2X2\]
Keterangan:
Mencari Eigen dengan R
# membuat data dummy
RNGkind(sample.kind = "Rounding")
set.seed(100)
x <- runif(200) # buat data random dengan distribusi normal
data <- data.frame(x=x, y=-x+runif(100, 1.05, 1.25)) # buat data frame
data <- scale(data) # agar skala antar variable sama
plot(data)# matrix covariance
A <- cov(data)
A#> x y
#> x 1.000000 -0.982719
#> y -0.982719 1.000000
# menghitung eigen
eigen(A)#> eigen() decomposition
#> $values
#> [1] 1.98271899 0.01728101
#>
#> $vectors
#> [,1] [,2]
#> [1,] -0.7071068 -0.7071068
#> [2,] 0.7071068 -0.7071068
Keterangan:
eigen(A)$values: Eigen value untuk tiap PC, besar variansi yang dapat ditangkap oleh tiap PC. Eigen value tertinggi adalah milik PC1, kedua tertinggi milik PC2, dan seterusnya.
eigen(A)$vectors: Eigen vector untuk tiap PC. Eigen vector untuk PC1 adalah kolom pertama [,1]; untuk PC2 adalah kolom kedua [,2].
Supervised vs Unsupervised
Dimensionality Reduction dengan PCA
Konsep:
Tujuan dimensionality reduction: mereduksi dimensi dengan tetap mempertahankan informasi sebesar mungkin
PCA merangkum variance (informasi) dari variabel awal dengan membuat sumbu baru (Principal Component), di mana informasi paling besar ditangkap oleh PC1, diikuti dengan PC2, dst
Banyaknya PC = Banyaknya variabel awal
PCA baik digunakan untuk data numerik yang saling berkorelasi. Dihasilkan nilai PC yang saling tidak berkorelasi
Istilah:
Eigen value = variance (informasi) yang ditangkap oleh tiap PC
Eigen vector = matrix rotasi untuk memproyeksikan data asli menjadi nilai PC
Workflow PCA
prcomp():pca$sdev adalah akar dari eigen valuepca$rotation adalah eigen vector tiap PCpca$x adalah nilai di tiap PC untuk tiap barissummary(prcomp())PCA tidak hanya berguna untuk dimensionality reduction namun baik untuk visualisasi high-dimensional data. Visualisasi dapat menggunakan biplot yang menampilkan:
Mari kita buat biplot dari 100 observasi pertama dari ppt
ppt_100 <- ppt %>%
head(100)
ppt_100_pca <- prcomp(ppt_100, scale = T)
biplot(ppt_100_pca, scale = F, cex = 0.6)Sumbu: - bawah dan kiri (merepresentasikan nilai pada observasi setelah PCA) - kanan dan atas (representasi dari kontribusi tiap var thd PC)
cor(ppt$TOTAL.UNITS, ppt$SALE.PRICE)#> [1] 0.1282114
ppt_100_pca$rotation#> PC1 PC2 PC3 PC4 PC5
#> RESIDENTIAL.UNITS 0.08651275 -0.69669893 0.056332773 -0.02580072 -0.07933713
#> COMMERCIAL.UNITS -0.36216824 0.11559840 -0.510870887 0.70840338 -0.05827666
#> TOTAL.UNITS -0.02536814 -0.69022011 -0.106422059 0.20084969 -0.10147273
#> LAND.SQUARE.FEET -0.51496295 -0.11878379 0.082731500 -0.26392120 0.65625154
#> GROSS.SQUARE.FEET -0.54594011 -0.08982028 -0.034704137 -0.04134819 0.12206617
#> YEAR.BUILT -0.19872684 0.03722296 0.846440178 0.47929655 -0.10056451
#> SALE.PRICE -0.50791375 0.03579098 0.001430911 -0.39505584 -0.72410861
#> PC6 PC7
#> RESIDENTIAL.UNITS 0.007257894 -0.7049397335
#> COMMERCIAL.UNITS 0.205575380 -0.2167702992
#> TOTAL.UNITS 0.070179991 0.6753254472
#> LAND.SQUARE.FEET 0.462100337 0.0013679038
#> GROSS.SQUARE.FEET -0.822231506 -0.0016927761
#> YEAR.BUILT 0.053123359 0.0007867168
#> SALE.PRICE 0.245636709 0.0008912568
ppt_100_pca$x[96,]#> PC1 PC2 PC3 PC4 PC5 PC6
#> -1.063228195 -6.008847310 0.312745012 -0.378234131 0.155742967 0.808400921
#> PC7
#> 0.000950834
ppt[c(51,32,96, 54,6), ]Biplot menggunakan PC1 dan PC2 secara default, karena PC tersebut yang merangkum informasi terbanyak dari data sehingga mampu mewakili data kita.
biplot(ppt_100_pca, scale = F, cex = 0.6) Notes: - Angka yang besar (pada sumbu) untuk representasi dari nilai observasi/titik2 setelah di PCA (bawah kiri) - Angka yang kecil2 (pada sumbu) untuk representasi kontribusi variabel (kanan atas) - Untuk antar sumbu apabila berhimpit dan searah (korelasi kuat positif) - Untuk antar sumbu apabila berhimpit dan berbalik arah (korelasi kuat negatif) - Untuk antar sumbu ~ 90 derajat (siku2) (korelasi ~0)
ppt_100_pca$x[96:100,]#> PC1 PC2 PC3 PC4 PC5 PC6
#> 96 -1.06322820 -6.0088473 0.31274501 -0.37823413 0.1557430 0.80840092
#> 97 0.78792455 0.6673594 -0.57993705 0.03398660 -0.1205905 -0.13921107
#> 98 0.01244112 0.8733915 0.06854775 0.08514277 0.2634052 -0.07600245
#> 99 0.72967334 0.9496234 -0.35840961 -0.33698321 0.2116062 -0.10873887
#> 100 0.30938895 -2.0009468 -0.06438969 0.51766854 -0.3724173 0.06100310
#> PC7
#> 96 0.0009508340
#> 97 -0.0001513919
#> 98 0.0009363043
#> 99 0.0004405085
#> 100 -0.0003498938
ppt_100_pca$rotation#> PC1 PC2 PC3 PC4 PC5
#> RESIDENTIAL.UNITS 0.08651275 -0.69669893 0.056332773 -0.02580072 -0.07933713
#> COMMERCIAL.UNITS -0.36216824 0.11559840 -0.510870887 0.70840338 -0.05827666
#> TOTAL.UNITS -0.02536814 -0.69022011 -0.106422059 0.20084969 -0.10147273
#> LAND.SQUARE.FEET -0.51496295 -0.11878379 0.082731500 -0.26392120 0.65625154
#> GROSS.SQUARE.FEET -0.54594011 -0.08982028 -0.034704137 -0.04134819 0.12206617
#> YEAR.BUILT -0.19872684 0.03722296 0.846440178 0.47929655 -0.10056451
#> SALE.PRICE -0.50791375 0.03579098 0.001430911 -0.39505584 -0.72410861
#> PC6 PC7
#> RESIDENTIAL.UNITS 0.007257894 -0.7049397335
#> COMMERCIAL.UNITS 0.205575380 -0.2167702992
#> TOTAL.UNITS 0.070179991 0.6753254472
#> LAND.SQUARE.FEET 0.462100337 0.0013679038
#> GROSS.SQUARE.FEET -0.822231506 -0.0016927761
#> YEAR.BUILT 0.053123359 0.0007867168
#> SALE.PRICE 0.245636709 0.0008912568
Dari biplot di atas, outlier terletak pada observasi 51,96,32
# outlier: 96 dan 32
# non-outlier: 26 dan 27
ppt_100[c(96, 32, 26, 27),]Observasi 96 dan 32 merupakan outlier dan ternyata karena nilai land square feet, total units, RU, dan GSF tinggi
biplot(ppt_100_pca, scale = F, cex = 0.6)Loading score dilihat dari jarak titik nol (pusat) ke ujung panah, secara horizontal (PC1) / vertikal (PC2). Semakin jauh panah, semakin banyak informasi yang dirangkum. Bila panah kearah negatif, maka bila nilai variabel awal tinggi, nilai di PC akan rendah, dan sebaliknya.
Berdasarkan penjelasan tersebut, maka:
Dari panah merah tersebut, kita tahu variable mana yang paling banyak berkontribusi untuk tiap PC. Namun kita kesulitan untuk mengurutkan kontribusinya, mari kita gunakan fungsi fviz_contrib() untuk melihat urutan kontribusi variabel ke tiap PC
library(factoextra)
fviz_contrib(X = ppt_100_pca, axes = 1, choice = "var")Catatan: garis putus-putus merah adalah batas nilai kontribusi yang diharapkan apabila PC menangkap informasi dari tiap variable dengan sama rata. Cara menghitungnya adalah 100%/jumlah variable awal. Pada kasus ini 100/7 = 14.28%. Bila melebihi batas tersebut, maka dapat dikatakan variable tersebut berkontribusi terhadap PC tersebut.
propertyBuatlah biplot dengan menggunakan 300 data pertama dari ppt:
# your code here
ppt_300 <- ppt %>% head(300)
ppt_300_pca <- prcomp(ppt_300, scale = T) # untuk membuat PCA
biplot(ppt_300_pca, #object PCA
cex = 0.6, # besar/kecil huruf
scale = F) # krn scr default di biplot dilakukan scaling lagi, agar bisa diinterpretasi maka scale = FMengidentifikasi outlier - yang jauh dari kumpulan datanya - data yang bergerombol di titik 0, artinya mereka ada di pusat datanya (nilainya mendekati rata2 dari var tsb)
_ Arah panah yang searah dan berhimpit artinya antar variabel berkorelasi kuat positif, sedangkan kalau bertolak belakang artinya korelasi negatif
Misal cek kolom RU dan TU
cor(ppt_300) # matrix korelasi#> RESIDENTIAL.UNITS COMMERCIAL.UNITS TOTAL.UNITS
#> RESIDENTIAL.UNITS 1.000000000 0.001875965 0.9216324
#> COMMERCIAL.UNITS 0.001875965 1.000000000 0.3897906
#> TOTAL.UNITS 0.921632390 0.389790601 1.0000000
#> LAND.SQUARE.FEET 0.235597352 0.173947032 0.2844254
#> GROSS.SQUARE.FEET 0.254876895 0.170640932 0.3009420
#> YEAR.BUILT 0.023298887 0.241512190 0.1150792
#> SALE.PRICE 0.198690043 0.210051040 0.2644687
#> LAND.SQUARE.FEET GROSS.SQUARE.FEET YEAR.BUILT SALE.PRICE
#> RESIDENTIAL.UNITS 0.2355974 0.2548769 0.02329889 0.1986900
#> COMMERCIAL.UNITS 0.1739470 0.1706409 0.24151219 0.2100510
#> TOTAL.UNITS 0.2844254 0.3009420 0.11507919 0.2644687
#> LAND.SQUARE.FEET 1.0000000 0.9139556 0.34307429 0.8646037
#> GROSS.SQUARE.FEET 0.9139556 1.0000000 0.32366540 0.9151348
#> YEAR.BUILT 0.3430743 0.3236654 1.00000000 0.3427099
#> SALE.PRICE 0.8646037 0.9151348 0.34270990 1.0000000
library(GGally)
ggcorr(ppt_300, label = T) # visualisasiMisal kita fokus ke var total units:
ppt_300_pca$rotation#> PC1 PC2 PC3 PC4 PC5
#> RESIDENTIAL.UNITS 0.2737605 0.63165943 -0.23947396 0.208551857 -0.022430850
#> COMMERCIAL.UNITS 0.1947503 0.09852218 0.79609954 -0.491509703 0.034089262
#> TOTAL.UNITS 0.3276678 0.61995107 0.08836838 0.001214791 -0.007536551
#> LAND.SQUARE.FEET 0.4871950 -0.23130179 -0.16547808 -0.089263044 0.719929484
#> GROSS.SQUARE.FEET 0.4962735 -0.22281835 -0.18709983 -0.115986829 -0.021481525
#> YEAR.BUILT 0.2480002 -0.19981608 0.47153449 0.822024473 -0.011075096
#> SALE.PRICE 0.4843664 -0.25275563 -0.12780713 -0.133333377 -0.692383856
#> PC6 PC7
#> RESIDENTIAL.UNITS -0.020215753 -0.65138932326
#> COMMERCIAL.UNITS 0.022131779 -0.27451326956
#> TOTAL.UNITS -0.009637761 0.70734314597
#> LAND.SQUARE.FEET -0.394323628 0.00016188608
#> GROSS.SQUARE.FEET 0.809409529 -0.00023013987
#> YEAR.BUILT 0.019435782 0.00007212157
#> SALE.PRICE -0.433585012 0.00006214021
Misal kita fokus ke observasi 118:
ppt_300_pca$x[118,]#> PC1 PC2 PC3 PC4 PC5
#> 16.8447196432 -8.7647259959 -4.1000035322 -1.5908537942 -1.9570796486
#> PC6 PC7
#> 0.9761761542 -0.0002879833
Nilai pada PC1 dan PC2 sudah sesuai dengan sumbu kiri dan bawah
Panah representasi dari variabel pada data awal
Insight: Outlier
Dari biplot, identifikasilah:
# optional: untuk memvalidasi insight
ppt_300 %>%
mutate(id = 1:300) %>%
arrange(-SALE.PRICE)Insight: Variable Contribution
GROSS.SQUARE.FEET, LAND.SQUARE.FEET, dan SALE.PRICERESIDENTIAL.UNITS dan TOTAL.UNITSCOMMERCIAL.UNITS dan YEAR.BUILTUntuk memvalidasi insight di atas, visualisasikan kontribusi variabel tiap PC menggunakan fviz_contrib() dari library factoextra:
# your code here
library(factoextra)
# PC1
fviz_contrib(X = ppt_300_pca,
choice = "var",
axes = 1)# PC2
fviz_contrib(X = ppt_300_pca,
choice = "var",
axes = 2)Insight: Variable Correlation
Gunakan biplot, identifikasilah:
GROSS.SQUARE.FEET, LAND.SQUARE.FEET, dan SALE.PRICERESIDENTIAL.UNITS dan TOTAL.UNITSSALE.PRICE dan TOTAL.UNITS# opsional: validasi korelasi
ggcorr(ppt_300, label = T)USArrestsLakukan PCA untuk built-in dataset USArrests dan visualisasikan menggunakan fancy_biplot()!
# get function & data
source("R/biplot.R")
data("USArrests")
head(USArrests)Dataset USArrests merupakan data tentang jumlah tiga tindak kejahatan (Murder, Assault, dan Rape) per 100 ribu penduduk yang berhasil dicatat di 50 negara bagian US pada tahun 1973.
Deskripsi data:
Murder: murder arrests per 100,000Assault: assault arrests per 100,000UrbanPop: percent of population in urban areasRape: rape arrests per 100,000# PCA & Visualization
fancy_biplot(prcomp(USArrests, scale = T))Insight dari biplot:
Perhatikan cluster di mana Arizona, Colorado, Illinois, dan Texas berada. Negara bagian mana lagi yang ada dalam cluster tersebut? New York, dimana cluster tersebut tinggi untuk tingkat Rape
Di antara Ohio, Louisiana, North Dakota, dan Colorado, manakah yang memiliki:
Murder tertinggi: LouisinaAssault tertinggi: LouisinaRape tertinggi: ColoradoBiplot
Mengetahui variabel apa saja yang dirangkum oleh tiap PC (variable contribution)
Korelasi antar variable, dilihat dari sudut antar panah:
Letak titik observasi terhadap arah panah:
FactoMineR packageFactoMineR adalah R package yang dibuat untuk exploratory multivariate data analysis. Mari kita lakukan analisis PCA untuk Dataset Loan dari sebuah platform di quarter ke-4 tahun 2017.
loan <- read.csv("data_input/loan2017Q4.csv", stringsAsFactors = T) # mengubah chr menjadi fct
glimpse(loan)#> Rows: 1,556
#> Columns: 16
#> $ initial_list_status <fct> w, f, w, w, w, w, w, w, w, w, w, w, w, f, w, w, w,~
#> $ purpose <fct> debt_consolidation, debt_consolidation, debt_conso~
#> $ int_rate <dbl> 14.08, 9.44, 28.72, 13.59, 15.05, 10.91, 15.05, 10~
#> $ installment <dbl> 675.99, 480.08, 1010.30, 484.19, 476.33, 130.79, 3~
#> $ annual_inc <dbl> 156700, 50000, 25000, 175000, 109992, 49000, 65000~
#> $ dti <dbl> 19.11, 19.35, 65.58, 12.60, 10.00, 5.12, 22.38, 33~
#> $ verification_status <fct> Source Verified, Not Verified, Verified, Not Verif~
#> $ grade <fct> C, B, F, C, C, B, C, B, D, D, F, C, C, E, B, C, C,~
#> $ revol_bal <int> 21936, 5457, 23453, 31740, 2284, 2016, 14330, 2758~
#> $ inq_last_12m <int> 3, 1, 0, 0, 3, 5, 0, 1, 8, 1, 0, 12, 4, 8, 1, 3, 0~
#> $ delinq_2yrs <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,~
#> $ home_ownership <fct> MORTGAGE, RENT, OWN, MORTGAGE, MORTGAGE, MORTGAGE,~
#> $ not_paid <int> 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1,~
#> $ log_inc <dbl> 11.962088, 10.819778, 10.126631, 12.072541, 11.608~
#> $ verified <int> 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,~
#> $ grdCtoA <int> 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,~
Data berisikan 1556 customer dengan 16 variabel, dengan kolom not_paid yang mengindikasikan apakah customer tersebut gagal bayar pinjaman atau tidak. Berikut adalah deskripsi data yang lengkap:
initial_list_status: Either w (whole) or f (fractional). This variable indicates if the loan was a whole loan or fractional loan. For background: Some institutional investors have a preference to purchase loans in their entirety to obtain legal and accounting treatment specific to their situation - with the added benefit of “instant funding” to borrowerspurpose: Simplified from the original data; One of: credit_card, debt_consolidation, home_improvement, major_purchase and small_businessint_rate: Interest rate in percentagesinstallment: Monthly payment owed by the borrowerannual_inc: Self-reported annual income provided by the borrower / co-borrowers during applicationdti: A ratio of the borrower’s total monthly debt payments on his/her total obligations to the self-reported monthly income (debt to income ratio)verification_status: is the reported income verified, not verified, or if the income source was verifiedgrade: software-assigned loan graderevol_bal: total credit revolving balance (in the case of credit card, it refers to the portion of credit card spending that goes unpaid at the end of a billing cycle)inq_last_12m: number of credit inquiries in the last 12 monthsdelinq_2yrs: number of 30+ days past-due incidences of delinquency in the borrower’s credit file for the past 2 yearshome_ownership: one of MORTGAGE, OWN and RENTnot_paid: 1 for charged-off, past-due / grace period or defaulted, 0 for fully-paid loanslog_inc: log of annual_incverified: 0 for “Not verified” under verification_status, 1 otherwisegrdCtoA: 0 for a grade of A, B or C, 1 otherwiseAdakah tipe data yang belum tepat?
loan_clean <- loan %>%
mutate_at(.vars = c("not_paid", "grdCtoA", "verified"), as.factor)FactoMineR menyediakan dua fungsi utama dalam PCA:
PCA() untuk analisis PCAplot.PCA() untuk visualisasi💡 Kelebihan FactoMineR dari fungsi base biplot() adalah:
individual factor map dan variables factor mapnames(loan_clean)#> [1] "initial_list_status" "purpose" "int_rate"
#> [4] "installment" "annual_inc" "dti"
#> [7] "verification_status" "grade" "revol_bal"
#> [10] "inq_last_12m" "delinq_2yrs" "home_ownership"
#> [13] "not_paid" "log_inc" "verified"
#> [16] "grdCtoA"
Cara 1
# indeks kolom numerik
quantivar <- c(3:6,9:11,14) # var numerik
# indeks kolom kategorik
qualivar <- c(1,2,7,8,12,13,15,16) # var categoricalSebagai alternatif, apabila tidak ingin manual:
# nama kolom numerik
quanti <- loan_clean %>%
select_if(is.numeric) %>%
colnames()
# indeks kolom numerik
quantivar <- which(colnames(loan) %in% quanti)
# nama kolom kategorik
quali <- loan_clean %>%
select_if(is.factor) %>%
colnames()
qualivar <- which(colnames(loan) %in% quali)Buat PCA dengan PCA()
Parameter:
X: dataframe awalscale.unit = T: untuk melakukan scalingquali.sup: index kolom variable kategorikal (qualitative)graph = F: untuk tidak menampilkan plot secara langsungncp: banyaknya PC yang digunakan, default = 5# PCA using FactoMineR
library(FactoMineR)
loan_pca <- PCA(X = loan_clean, scale.unit = T, quali.sup = qualivar, graph = F, ncp = 8)Notes: ncp = 8 karena jumlah variabel numeric pada loan_clean adalah 8 (jumlah quantivar), tujuannya adalah kita bisa akses seluruh hasil PC
# cek nilai di tiap PC, ekuivalen dengan obj_pca$x
head(loan_pca$ind$coord)#> Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
#> 1 1.8726407 -0.02630256 -0.2620212 -0.3818516 -0.1440523 0.29379806
#> 2 -0.7953612 -0.59637757 0.1002547 0.9009347 0.4275856 0.77186601
#> 3 -1.3593014 4.73853507 -0.5401945 0.8218723 -0.5799966 0.74642058
#> 4 2.0152497 -0.64397562 -0.9941929 0.4797471 -0.3253111 -0.52434576
#> 5 0.6545361 -0.71599161 0.3267262 -0.5351297 -0.6294665 -0.01941924
#> 6 -0.9730541 -1.37420746 0.4126093 -1.2286800 0.1841724 -0.34258151
#> Dim.7 Dim.8
#> 1 0.4142317 0.2487058
#> 2 -0.1167101 -0.1040927
#> 3 0.8142216 -0.2268024
#> 4 0.5452035 0.2257019
#> 5 0.3436425 0.2451064
#> 6 -0.5383739 -0.1755210
Biplot dapat divisualisasikan menggunakan plot.PCA()
Plot sebaran observasi untuk mengetahui index yang dianggap sebagai outlier
Parameter:
x: object PCA (FactoMineR)choix = "ind": plot individual factor mapinvisible = "quali": menghilangkan label kolom kategorik, karena mengganggu visualselect = "contrib n": menampilkan indeks dari \(n\) outlier terluarhabillage: mewarnai titik berdasarkan index kolom/ nama kolom (bisa kategorikal maupun numerik)plot.PCA(x = loan_pca, # obj PCA
choix = "ind", # pilihan, kalau mau liat per observasi pakai ind
select = "contrib 5", # untuk identifikasi 5 outlier terluar
habillage = "purpose", # untuk mewarnai titik observasi (ditulis index/nama kolom)
invisible = "quali") # untuk menghilangkan garis panah cek data-data yang termasuk outlier
loan_clean[c(1228,1146,368,749,512),]cek output dr PCA (sama dengan obj_pca$x pada base R)
loan_pca$ind$coord[749,]#> Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 Dim.7
#> 12.9274431 -3.0905871 -0.8027414 -0.7007887 -0.2127467 -1.2044792 8.6655309
#> Dim.8
#> -7.3282417
Lima outlier terluar: 1228,1146,368,749,512
Parameter:
x: object PCA (FactoMineR)choix = "var": plot variable factor mapplot.PCA(x = loan_pca, choix = "var")utk cek nilai eigen vectors (besarnya kontribusi thd PC tertentu)
loan_pca$var$coord#> Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
#> int_rate 0.03654542 0.68194854 0.43163633 0.02695565 -0.36883724
#> installment 0.56432369 0.43476671 0.11703206 0.15624487 -0.34570189
#> annual_inc 0.88930271 -0.12014772 -0.05726255 -0.05053072 0.06550692
#> dti -0.25557584 0.73190607 -0.27580312 -0.03837392 0.38174826
#> revol_bal 0.55136181 0.37094415 -0.42453389 0.14683011 0.32889888
#> inq_last_12m 0.18305771 0.10857353 0.49152494 -0.75816324 0.32939466
#> delinq_2yrs 0.05152860 -0.06936228 0.63030805 0.62182210 0.45475297
#> log_inc 0.91272785 -0.21604640 0.01121066 -0.03093873 -0.01884030
#> Dim.6 Dim.7 Dim.8
#> int_rate -0.44257783 0.1195321468 0.019238287
#> installment 0.53467400 -0.2170534397 -0.043697024
#> annual_inc -0.06727446 0.3358663756 -0.259325342
#> dti 0.23656179 0.3407151270 0.060618128
#> revol_bal -0.31469567 -0.3864379359 -0.008388692
#> inq_last_12m 0.06069730 -0.1611694076 -0.011524710
#> delinq_2yrs 0.04201032 0.0005436539 -0.004156749
#> log_inc -0.00551407 0.1633062154 0.303503982
Insight:
cek korelasi pada data asli
cor(loan_clean %>% select_if(is.numeric))#> int_rate installment annual_inc dti revol_bal
#> int_rate 1.00000000 0.23592630 -0.03474266 0.16609373 0.06544226
#> installment 0.23592630 1.00000000 0.31483666 0.05361748 0.24796176
#> annual_inc -0.03474266 0.31483666 1.00000000 -0.18968136 0.37774994
#> dti 0.16609373 0.05361748 -0.18968136 1.00000000 0.16097269
#> revol_bal 0.06544226 0.24796176 0.37774994 0.16097269 1.00000000
#> inq_last_12m 0.10461178 0.04363990 0.12626494 0.01070311 -0.02716944
#> delinq_2yrs 0.05706945 0.03516166 0.01486787 -0.07816637 -0.03743149
#> log_inc -0.07522290 0.37247865 0.81384981 -0.32775948 0.34368507
#> inq_last_12m delinq_2yrs log_inc
#> int_rate 0.104611775 0.057069448 -0.0752229
#> installment 0.043639902 0.035161663 0.3724787
#> annual_inc 0.126264937 0.014867866 0.8138498
#> dti 0.010703114 -0.078166369 -0.3277595
#> revol_bal -0.027169436 -0.037431491 0.3436851
#> inq_last_12m 1.000000000 -0.007425338 0.1362335
#> delinq_2yrs -0.007425338 1.000000000 0.0398727
#> log_inc 0.136233530 0.039872705 1.0000000
Untuk lebih jelas dan objektif, kontribusi/korelasi tiap variable ke tiap PC dapat dilihat menggunakan dimdesc(). Semakin tinggi nilai korelasi, semakin banyak informasi yang dirangkum pada PC tersebut.
# dimdisc: dimension description
dim <- dimdesc(loan_pca)
# variable yang berkontribusi untuk PC1
as.data.frame(dim$Dim.1$quanti)💡 Mirip dengan fviz_contrib() yang kita pelajari sebelumnya
library(factoextra)
fviz_contrib(loan_pca, choice = "var")⚠️ Besar kontribusi tidak dilihat dari tanda pada korelasi, melainkan nilai absolutnya
Menampilkan cumulative proportion dari PCA:
# PCA summary
summary(loan_pca) # menampilkan summary keseluruhan#>
#> Call:
#> PCA(X = loan_clean, scale.unit = T, ncp = 8, quali.sup = qualivar,
#> graph = F)
#>
#>
#> Eigenvalues
#> Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 Dim.7
#> Variance 2.349 1.405 1.099 1.013 0.829 0.647 0.492
#> % of var. 29.365 17.563 13.732 12.664 10.368 8.084 6.153
#> Cumulative % of var. 29.365 46.929 60.661 73.325 83.693 91.777 97.931
#> Dim.8
#> Variance 0.166
#> % of var. 2.069
#> Cumulative % of var. 100.000
#>
#> Individuals (the 10 first)
#> Dist Dim.1 ctr cos2 Dim.2 ctr cos2
#> 1 | 2.016 | 1.873 0.096 0.863 | -0.026 0.000 0.000 |
#> 2 | 1.616 | -0.795 0.017 0.242 | -0.596 0.016 0.136 |
#> 3 | 5.184 | -1.359 0.051 0.069 | 4.739 1.027 0.835 |
#> 4 | 2.534 | 2.015 0.111 0.632 | -0.644 0.019 0.065 |
#> 5 | 1.382 | 0.655 0.012 0.224 | -0.716 0.023 0.269 |
#> 6 | 2.233 | -0.973 0.026 0.190 | -1.374 0.086 0.379 |
#> 7 | 1.171 | -0.535 0.008 0.209 | -0.001 0.000 0.000 |
#> 8 | 2.232 | 0.934 0.024 0.175 | 0.871 0.035 0.152 |
#> 9 | 3.043 | 2.215 0.134 0.530 | 0.371 0.006 0.015 |
#> 10 | 1.004 | 0.076 0.000 0.006 | 0.053 0.000 0.003 |
#> Dim.3 ctr cos2
#> 1 -0.262 0.004 0.017 |
#> 2 0.100 0.001 0.004 |
#> 3 -0.540 0.017 0.011 |
#> 4 -0.994 0.058 0.154 |
#> 5 0.327 0.006 0.056 |
#> 6 0.413 0.010 0.034 |
#> 7 -0.738 0.032 0.397 |
#> 8 -1.123 0.074 0.253 |
#> 9 0.903 0.048 0.088 |
#> 10 -0.174 0.002 0.030 |
#>
#> Variables
#> Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3
#> int_rate | 0.037 0.057 0.001 | 0.682 33.098 0.465 | 0.432
#> installment | 0.564 13.556 0.318 | 0.435 13.453 0.189 | 0.117
#> annual_inc | 0.889 33.665 0.791 | -0.120 1.027 0.014 | -0.057
#> dti | -0.256 2.780 0.065 | 0.732 38.125 0.536 | -0.276
#> revol_bal | 0.551 12.941 0.304 | 0.371 9.793 0.138 | -0.425
#> inq_last_12m | 0.183 1.426 0.034 | 0.109 0.839 0.012 | 0.492
#> delinq_2yrs | 0.052 0.113 0.003 | -0.069 0.342 0.005 | 0.630
#> log_inc | 0.913 35.462 0.833 | -0.216 3.322 0.047 | 0.011
#> ctr cos2
#> int_rate 16.959 0.186 |
#> installment 1.247 0.014 |
#> annual_inc 0.298 0.003 |
#> dti 6.924 0.076 |
#> revol_bal 16.405 0.180 |
#> inq_last_12m 21.991 0.242 |
#> delinq_2yrs 36.163 0.397 |
#> log_inc 0.011 0.000 |
#>
#> Supplementary categories (the 10 first)
#> Dist Dim.1 cos2 v.test Dim.2 cos2 v.test
#> f | 0.279 | -0.098 0.123 -1.292 | 0.173 0.382 2.939 |
#> w | 0.073 | 0.026 0.123 1.292 | -0.045 0.382 -2.939 |
#> credit_card | 0.330 | -0.064 0.037 -0.836 | -0.189 0.331 -3.213 |
#> debt_consolidation | 0.141 | -0.019 0.018 -0.645 | 0.120 0.719 5.319 |
#> home_improvement | 0.445 | 0.107 0.057 0.932 | -0.276 0.384 -3.117 |
#> major_purchase | 0.807 | 0.405 0.252 1.918 | -0.274 0.115 -1.678 |
#> small_business | 0.358 | 0.051 0.020 0.173 | -0.008 0.000 -0.034 |
#> Not Verified | 0.385 | -0.176 0.209 -3.476 | -0.287 0.555 -7.328 |
#> Source Verified | 0.184 | 0.094 0.259 1.749 | -0.113 0.380 -2.737 |
#> Verified | 0.572 | 0.117 0.042 1.882 | 0.515 0.811 10.753 |
#> Dim.3 cos2 v.test
#> f 0.140 0.252 2.699 |
#> w -0.037 0.252 -2.699 |
#> credit_card -0.234 0.505 -4.493 |
#> debt_consolidation 0.030 0.044 1.496 |
#> home_improvement 0.137 0.095 1.752 |
#> major_purchase 0.321 0.158 2.224 |
#> small_business 0.261 0.533 1.305 |
#> Not Verified -0.151 0.153 -4.358 |
#> Source Verified 0.054 0.087 1.478 |
#> Verified 0.132 0.053 3.116 |
loan_pca$eig # menampilkan cummulative proportion variance/eigen value#> eigenvalue percentage of variance cumulative percentage of variance
#> comp 1 2.3492124 29.365155 29.36515
#> comp 2 1.4050728 17.563410 46.92857
#> comp 3 1.0985925 13.732406 60.66097
#> comp 4 1.0131555 12.664444 73.32542
#> comp 5 0.8294041 10.367552 83.69297
#> comp 6 0.6467516 8.084395 91.77736
#> comp 7 0.4922722 6.153403 97.93076
#> comp 8 0.1655389 2.069236 100.00000
ini kalau di base R (untuk cek proporsi variancenya)
summary(pca)#> Importance of components:
#> PC1 PC2 PC3 PC4 PC5 PC6 PC7
#> Standard deviation 1.6918 1.0941 1.0073 0.9995 0.8426 0.46601 0.005136
#> Proportion of Variance 0.4089 0.1710 0.1449 0.1427 0.1014 0.03102 0.000000
#> Cumulative Proportion 0.4089 0.5799 0.7248 0.8675 0.9690 1.00000 1.000000
Misalkan kita ingin mempertahankan informasi sebesar 80%, maka berapa PC yang kita gunakan? PC 1-5
# mengambil data hasil PCA sebanyak PC yang dibutuhkan
loan_keep <- loan_pca$ind$coord[,1:5] %>%
as.data.frame()loan_keep kemudian dapat dikembalikan ke data frame awal (menggantikan variable numerik awal) dan digunakan untuk supervised learning.
Notes: Umumnya yang digunakan adalah 70-90% informasi
Fungsi reconst() memungkinkan kita untuk melakukan rekonstruksi data hasil PCA ke data semula. Umumnya rekonstruksi ini digunakan pada data gambar untuk kasus image compression.
Misalnya gambar asli adalah 40x40 pixel (1600 kolom), lalu direduksi dengan menggunakan 50 PC saja. Untuk menampilkan data berupa gambar, tetap diperlukan struktur data yang 1600 kolom yang mewakili tiap pixel. Kita bisa melakukan rekonstruksi gambar dari 50 PC ke 1600 kolom kembali. Namun, nilai hasil rekonstruksi tersebut menjadi sedikit berbeda dari nilai awal.
Catatan: rekonstruksi data jarang dilakukan pada data tabular, karena nilai hasil rekonstruksi tidak merepresentasikan nilai asli, sehingga kurang ada manfaatnya.
Berikut contoh rekonstruksi data dengan menggunakan data loan:
# reconstruct data menggunakan PC1 - PC5
loan_reconst <- reconst(loan_pca,
ncp = 5) # jumlah PC yang ingin di reconstruct
head(loan_reconst, 3)#> int_rate installment annual_inc dti revol_bal inq_last_12m delinq_2yrs
#> 1 14.53469 663.7699 155688.191 14.89860 29305.69 3.2019752 -0.01035623
#> 2 12.01608 310.8581 53601.426 17.34098 10443.12 0.7778570 0.96742349
#> 3 30.34033 929.1623 -6418.011 58.04061 39025.64 0.3136626 -0.03308597
#> log_inc
#> 1 11.79575
#> 2 10.88521
#> 3 10.11759
Bandingkan dengan nilai asli sebelum dilakukan PCA:
# data awal (coba bandingkan dengan hasil reconstruct di atas)
loan_clean %>%
select_if(is.numeric) %>%
head(3)Kita juga dapat rekonstruksi data hingga 100% sama persis bila menggunakan seluruh PC, dengan menggunakan parameter ncp = 8 di fungsi reconst().
Bagian berikut menyajikan contoh rekonstruksi PCA pada data gambar.
Pada kasus face recognition, data gambar yang memiliki banyak sekali dimensi dapat mengalami apa yang disebut sebagai the curse of dimensionality.
Setelah melewati batas tertentu, kemampuan suatu model dalam memprediksi hasil akan semakin menurun seiring dengan bertambahnya dimensi.
Dengan melakukan dimensionality reduction (PCA) pada data gambar dan merekonstruksi data gambar dari data hasil PCA, kita dapat meningkatkan performa model face recognition kita.
# read `faceData` 32 pixel
load("data_input/face.rda")
# PCA
face_pca <- PCA(faceData, graph = F, ncp = 32)
head(face_pca$eig, 10)#> eigenvalue percentage of variance cumulative percentage of variance
#> comp 1 12.6194690 39.4358407 39.43584
#> comp 2 7.3922489 23.1007778 62.53662
#> comp 3 4.7536745 14.8552327 77.39185
#> comp 4 2.2909120 7.1590999 84.55095
#> comp 5 1.2400024 3.8750076 88.42596
#> comp 6 0.7862118 2.4569118 90.88287
#> comp 7 0.6251630 1.9536343 92.83650
#> comp 8 0.5077062 1.5865819 94.42309
#> comp 9 0.4519846 1.4124518 95.83554
#> comp 10 0.2936547 0.9176711 96.75321
# create reconstructed data
face_recon <- reconst(face_pca, ncp = 2)
face_recon2 <- reconst(face_pca, ncp = 4)
face_recon3 <- reconst(face_pca, ncp = 6)
face_recon4 <- reconst(face_pca, ncp = 32) # mengambil semua PC# make function to visualize image data
showMatrix <- function(x, title){
image(t(x[nrow(x):1,]),
xaxt = 'n', yaxt = 'n',
col = gray((0:32)/32),
main = title,
font.main=4,
cex.main=0.5
)
}
# visualize image data
par(mfrow=c(2,2), mar=c(0.5,0.5,1.5,0.5))
showMatrix(faceData, title = 'Original Image')
showMatrix(face_recon, title = 'Reconstructed: 2 dimensions')
showMatrix(face_recon3, title = 'Reconstructed: 6 dimensions')
showMatrix(face_recon4, title = 'Reconstructed: All PCs')Fungsi-fungsi
biplot(prcomp()): membuat biplot dengan base RPCA(): membuat objek PCA dari package FactorMineRplot.PCA(): membuat biplot dari object PCA():
choix = "ind": individual factor map, terdapat habillage untuk mewarnai observasi berdasarkan index kolom,choix = "var": variables factor mapdimdesc(): dimension description, untuk melihat korelasi (kontribusi) variable ke tiap PCreconst(): rekonstruksi data, umumnya untuk data gambarClustering adalah pengelompokan data berdasarkan karakteristiknya. Clustering bertujuan untuk menghasilkan cluster dimana:
K-means adalah salah satu algoritma clustering yang centroid-based, artinya tiap cluster memiliki satu centroid yang mewakili cluster tersebut.
💡 Banyaknya cluster \(k\) ditentukan oleh user.
# pakai 'interval' yang lebih tinggi bila animasi terlalu cepat
# jalankan command ini di console:
library(animation)
RNGkind(sample.kind = "Rounding")
set.seed(100)
ani.options(interval = 2)
par(mar = c(3, 3, 1, 1.5), mgp = c(1.5, 0.5, 0))
kmeans.ani()💡 Kunjungi link berikut untuk visualisasi yang interaktif: https://www.naftaliharris.com/blog/visualizing-k-means-clustering/
Kita sebagai seorang data scientist sebuah toko whisky diminta untuk membuat product recommendation untuk whisky berdasarkan preferensi rasa masing-masing customer!
Tujuan: membentuk kelompok whisky yang memiliki karakteristik rasa khas pada tiap clusternya
Data yang digunakan berupa data penyulingan Malt Whisky dari 86 pabrik penyulingan, diperoleh dari penelitian Dr. Wisehart (Universitas St. Andrews). Setiap whisky diberi skor 0-4 dari 12 kategori cita rasa berdasarkan uji organoleptik:
whisky <- read.delim("data_input/whiskies.txt", sep = ",")
head(whisky)library(tidyverse)
# meng-assign nilai dari kolom Distillery menjadi rownames
whisky_clean <- whisky %>%
column_to_rownames(var = "Distillery")
# membuang kolom yang tidak digunakan
whisky_clean <- whisky_clean %>%
select(-c(RowID, Postcode, Latitude, Longitude))
head(whisky_clean)Cek missing values:
anyNA(whisky_clean)#> [1] FALSE
Cek skala antar variabel
summary(whisky_clean)#> Body Sweetness Smoky Medicinal
#> Min. :0.00 Min. :1.000 Min. :0.000 Min. :0.0000
#> 1st Qu.:2.00 1st Qu.:2.000 1st Qu.:1.000 1st Qu.:0.0000
#> Median :2.00 Median :2.000 Median :1.000 Median :0.0000
#> Mean :2.07 Mean :2.291 Mean :1.535 Mean :0.5465
#> 3rd Qu.:2.00 3rd Qu.:3.000 3rd Qu.:2.000 3rd Qu.:1.0000
#> Max. :4.00 Max. :4.000 Max. :4.000 Max. :4.0000
#> Tobacco Honey Spicy Winey
#> Min. :0.0000 Min. :0.000 Min. :0.000 Min. :0.0000
#> 1st Qu.:0.0000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:0.0000
#> Median :0.0000 Median :1.000 Median :1.000 Median :1.0000
#> Mean :0.1163 Mean :1.244 Mean :1.384 Mean :0.9767
#> 3rd Qu.:0.0000 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:1.0000
#> Max. :1.0000 Max. :4.000 Max. :3.000 Max. :4.0000
#> Nutty Malty Fruity Floral
#> Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
#> 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000
#> Median :2.000 Median :2.000 Median :2.000 Median :2.000
#> Mean :1.465 Mean :1.802 Mean :1.802 Mean :1.698
#> 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000
#> Max. :4.000 Max. :3.000 Max. :3.000 Max. :4.000
Diskusi: Pada data whisky, apakah skala nilai antar variable berbeda? Apakah perlu dilakukan scaling?
Tidak perlu, karena di deskripsi datanya, untuk tiap rasa diberi skala yang sama, yakni 0-4
Parameter:
x: datasetcenters: banyaknya centroid \(k\)Note: perlu dilakukan set.seed() karena terdapat random initialization pada tahap awal k-means
# k-means dengan 3 cluster
RNGkind(sample.kind = "Rounding")
set.seed(100)
whi_km <- kmeans(x = whisky_clean, # data yang sudah dicleaning
centers = 3) # jumlah k (cluster yang diinginkan)Hasil dari k-means:
whi_km$iter#> [1] 2
whi_km$size#> [1] 34 11 41
whi_km$centers#> Body Sweetness Smoky Medicinal Tobacco Honey Spicy Winey
#> 1 2.500000 2.323529 1.588235 0.1764706 0.05882353 1.8823529 1.647059 1.6764706
#> 2 2.909091 1.545455 2.909091 2.7272727 0.45454545 0.4545455 1.454545 0.5454545
#> 3 1.487805 2.463415 1.121951 0.2682927 0.07317073 0.9268293 1.146341 0.5121951
#> Nutty Malty Fruity Floral
#> 1 1.823529 2.088235 1.911765 1.7058824
#> 2 1.545455 1.454545 1.181818 0.5454545
#> 3 1.146341 1.658537 1.878049 2.0000000
head(whi_km$cluster) # menghasilkan cluster untuk tiap rownames (nama-nama whiskynya)#> Aberfeldy Aberlour AnCnoc Ardbeg Ardmore ArranIsleOf
#> 1 1 3 2 1 3
Kebaikan hasil clustering dapat dilihat dari 3 nilai:
Within Sum of Squares ($withinss): jumlah jarak kuadrat dari tiap observasi ke centroid tiap cluster.
Between Sum of Squares ($betweenss): jumlah jarak kuadrat terbobot dari tiap centroid ke rata-rata global. Dibobotkan berdasarkan banyaknya observasi pada cluster.
Total Sum of Squares ($totss): jumlah jarak kuadrat dari tiap observasi ke rata-rata global.
Ilustrasi lihat pada slide: https://docs.google.com/presentation/d/1XzKN4fYqQ1VnS9x9A3GvkVi0torLmpfIGXwVeYzMkOc/edit#slide=id.gb7882dd136_0_1472
# nilai WSS dan BSS/TSS
whi_km$withinss#> [1] 179.20588 65.45455 202.53659
whi_km$betweenss#> [1] 218.6402
whi_km$tot.withinss#> [1] 447.197
Clustering yang “baik”:
WSS semakin rendah: jarak observasi di 1 kelompok yang sama semakin rendah, artinya tiap cluster memiliki karakteristik yang semakin mirip (WSS ~ 0)
Tujuan: WSS semakin kecil dan Jumlah cluster juga kecil
\(\frac{BSS}{TSS}\) mendekati 1, karena kelompok hasil clustering semakin mewakili persebaran data yang sesungguhnya
whi_km$betweenss/whi_km$totss#> [1] 0.3283688
Diskusi: Coba buatlah k-means dengan \(k\) sebesar mungkin, misalnya 80. Kemudian cek WSS dan BSS/TSSnya, apakah clustering yang terbentuk dapat dikatakan ideal?
RNGkind(sample.kind = "Rounding")
set.seed(123)
# clustering
whi_km_80 <- kmeans(whisky_clean, centers = 80)
# nilai WSS dan BSS/TSS
whi_km_80$withinss#> [1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0
#> [20] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0
#> [39] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.5 0.0 0.0 0.0 0.0
#> [58] 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> [77] 0.0 0.0 0.0 0.0
whi_km_80$betweenss/whi_km_80$totss#> [1] 0.9917397
Jawaban: …
Semakin tinggi \(k\):
Kalau begitu apakah kita selalu memilih \(k\) = banyak observasi? Bagaimana menentukan \(k\) optimum?
fviz_nbclust()library(factoextra)
fviz_nbclust(x = whisky_clean, FUNcluster = kmeans, method = "wss")Pilih nilai \(k\) dimana ketika \(k\) ditambah, penurunan Total WSS tidak terlalu drastis (atau dapat dikatakan melandai).
Pilih k = 5 (karena dari cluster 5 ke 6 penurunan nilai total withinss landai/tidak signifikan)
Buat ulang k-means dengan \(k\) optimum:
RNGkind(sample.kind = "Rounding")
set.seed(100)
whi_km_opt <- kmeans(whisky_clean, centers = 5)# memasukkan label cluster ke data awal
whisky_clean$cluster <- as.factor(whi_km_opt$cluster)
head(whisky_clean)# melakukan profiling dengan summarise data
whisky_profile <- whisky_clean %>%
group_by(cluster) %>%
summarise_all(.funs = "mean")Alternatifnya, cukup lihat centroid: whi_km_opt$centers
Profiling:
# optional: mempermudah profiling
whisky_profile %>%
tidyr::pivot_longer(-cluster) %>%
group_by(name) %>%
summarize(cluster_min_val = which.min(value),
cluster_max_val = which.max(value))Misal ada 1 pelanggan yang menyukai whisky Laphroig datang ke toko kita, namun stok whisky tersebut sedang kosong. Kira-kira whisky apa yang akan kita rekomendasikan?
whisky_clean %>%
filter(rownames(whisky_clean) == "Laphroig")whisky_clean %>%
filter(cluster == 2)Misalkan ada kebutuhan untuk memvisualisasikan hasil clustering tersebut. Namun pada dataset whisky terdapat 12 variabel numerik, bagaimana cara kita memvisualisasikannya pada plot 2 dimensi?
# package factoextra
fviz_cluster(object = whi_km_opt,
data = whisky_clean %>% select(-cluster))Misal kita ingin buat sebuah visualisasi yang mempermudah cluster profiling, dimana tampil gabungan individual + variable factor map menjadi satu. Visualisasi dapat dibuat menggunakan fungsi fviz_pca_biplot() dari package factoextra
# PCA menggunakan FactoMineR
whisky_pca <- PCA(whisky_clean,
scale.unit = F, # tidak perlu scaling
quali.sup = 13, # kolom cluster
graph = F)# visualisasi PCA + hasil kmeans clustering
fviz_pca_biplot(whisky_pca,
habillage = 13, # kolom cluster
geom.ind = "point",
addEllipses = T)Clustering
Mengelompokkan data berdasarkan kemiripan karakteristik, harapannya:
K-means Clustering
Workflow K-means
kmeans(data, center = k)Inclass Questions:
http://www.gastonsanchez.com/visually-enforced/how-to/2012/10/13/MCA-in-R/